Practical solutions for perturbed f{R) gravity 
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We examine the accuracy of the quasi-static approximation and a parametric form commonly 
employed for evolving hnear perturbations in f{R) gravity theories and placing cosmological con- 
straints. In particular, we analyze the nature and the importance of the near horizon effects that 
are often neglected. We find that for viable models, with a small present value of the scalaron 
field, such corrections are entirely negligible with no impact on observables. We also find that the 
one-parameter form, commonly used to represent the modified equations for linear perturbations in 
f{R), leads to theoretical systematic errors that are relatively small and, therefore, is adequate for 
placing constraint on f{R) models. 



I. INTRODUCTION 

Understanding the physics behind the observed accel- 
erating expansion of the universe [1, 2] is one of the most 
important task of modern science. The simplest theoret- 
ical explanation within the context of General Relativity 
(GR) [3] is the cosmological constant A, which also hap- 
pens to be in good agreement with the current cosmic 
microwave background (CMB) [4], supenovae (SN) [5, 6], 
and baryon acoustic oscillations (BAO) [7, 8] data. How- 
ever, the observed A is many orders of magnitude smaller 
than the vacuum energy density predicted by quantum 
field theory [9], implying either an extraordinary de- 
gree of fine-tuning or an anthropic argument [10, 11]. 
Many alternatives to A have been proposed in the liter- 
ature so far, and they can be broadly divided into two 
classes: dynamical dark energy [12] and modified gravity 
(MG) [13, 14]. 

A very popular model of MG involves adding a function 
of the Ricci scalar, f{R), to the Einstein-Hilbert action 
of GR, such that the resulting equations for the metric 
admit self-accelerating solutions [15-19]. In fact, given 
the 4th order nature of the f{R) equations of motion, it is 
possible to design /(i?) functions to match any expansion 
history of the universe [20]. Prior to its appearance in 
the context of the current cosmic acceleration, f{R) was 
introduced by Starobinsky as the first working model of 
Inflation [21, 22]. 

The 4th order nature of the equations implies the exis- 
tence of an additional scalar degree of freedom coupled to 
all matter with gravitational strength^. The scalar field, 
dubbed "scalaron" , mediates an attractive force, with the 
trace of the energy-momentum acting as charge. It can 
evade solar system tests [25] constraining fifth force in- 
teractions thanks to the Chameleon mechanism [26] in 



^ By transforming to a conformally related frame f{R) can be 
explicitly written as a scalar-tensor gravity theory [23, 24]. 



which the scalaron is screened in high density environ- 
ments. While most choices of the f{R) function result 
in a time-varying equation of state w for the effective 
dark energy fluid, viable models are forced to have an 
expansion history that is practically indistinguishable 
from that of the A Cold Dark Matter (LCDM) model 
with w — —\ [17-19, 27-29]. There can, however, be de- 
tectable differences between /(i?) and LCDM when it 
comes to the dynamics of clustering of matter, both on 
linear and non-linear scales [17, 20, 30-37]. Searching for 
evidence of modified growth dynamics is indeed one of 
the main science goals of upcoming and future large scale 
surveys such as DES [38], Euclid [39] and LSST [40]. 

The dynamics of linear cosmological perturbations in 
f{R) has been studied in [20, 30, 41]. In principle, the 
non-linearity of the f{R) function means that it is not 
guaranteed that thepredictions of perturbation theory 
will necessarily agree with the actual dynamics of clus- 
tering on large scales. In fact, if all of the matter was 
in sufficiently dense form, the scalaron force would al- 
ways be suppressed by the Chameleon mechanism and 
the growth on linear scales would be identical to that in 
GR. In reality, most of the matter in the universe resides 
in diffuse regions that do not screen the fifth force, and 
N-body simulations [31, 32, 35, 42, 43] have confirmed 
that predictions of linear theory are indeed recovered on 
large scales. Given the degree of fine-tuning required to 
design viable f{R) models, they can hardly be considered 
a serious alternative to A. However, they provide a simple 
working toy-model for studying strongly coupled scalar- 
tensor theories with Chameleon screening that help us 
better understand ways in which gravity can be tested 
with upcoming observational data. 

Solving exact perturbed equations in MG models for 
the purpose of calculating observables can be challenging. 
Often authors rely on approximate parametrized forms of 
the equations which are meant to capture the main phys- 
ical features of the model while providing a simplified 
framework to evolve perturbations. A fairly general five- 
parameter model was introduced by Bertschinger and 
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Zukin (BZ) in [44]. A single parameter version of the 
BZ parametrization was used in [45-47] and other hter- 
ature to mimic the solutions of the perturbation equa- 
tions in f{R) theories when deriving cosmological con- 
straints. However, the accuracy of this single parameter 
form has not been consistently tested. In this work, we 
compare the single-parameter form used in [45-47] to the 
exact numerical solutions in f{R) and estimate the sys- 
tematic theoretical error. We analyze several aspects of 
this parametrization, including the accuracy of the quasi- 
static approximation on which it is based. We find that 
the near-horizon effects are more prominent for large val- 
ues of the scalaron today, which are already ruled, but 
are negligible for f{R) models in the viable range. Over- 
all, we conclude that the single parameter BZ form is an 
adequate representation of solutions in linearized f{R) 
models for the purpose of testing them with the most 
advanced future surveys. 



II. /(i?) THEORIES 

We consider f{R) models of gravity in the Jordan 
frame, for which the action reads: 



curvature regime and a non-tachyonic scalaron [27], and 
1 -|- //{ > to prevent the scalaron from turning into a 
ghost [48-50] . Considering that we want to reproduce GR 
at early times, i.e. we want — >■ as i? — oo, these con- 
ditions imply that fn will be a monotonically increasing 
function of R that asymptotes to zero from below. If one 
furthermore wants to satisfy existing constraints from lo- 
cal tests of gravity, has to be small at recent epochs; 
solar system and nearby universe tests impose a bound 
of l/fll < 10-6 [17, 51]. 

The Friedmann equation (2) can be re-written as a 
second order equation for /(a). Given an expansion his- 
tory H{a), and hence R{a), one can use (2) to solve for 
/(a) and then determine f{R) [20, 30]. In this "designer" 
approach, the amplitude of the decaying mode is set to 
zero, while the cosmologically relevant growing mode so- 
lution is not unique - for each expansion history there 
is a family of f{R) models parametrized by a boundary 
condition, such as the value of df /dR today, which we 
denote as In the literature, constraints on f{R) are 
often presented in terms of a parameter Bq [20] , which is 
related to the mass of the scalaron today. It is another 
way of setting the boundary condition that labels the 
models and is given by the value of the function 
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d*x^-g [R + f{R)] + S^[x^,g^.] , (1) 



where f{R) is a general function of the Ricci scalar. Mat- 
ter is minimally coupled and, therefore, the matter fields, 
Xi, fall along geodesies of the metric g^^. The field equa- 
tions obtained from varying action (1) with respect to 
(7^1, are of higher order than those in GR. At the back- 
ground level, the cosmological evolution is described by 
the Friedmann equation: 



o a 
and the trace equation: 



Ofn=l(R + 2f-Rf,)-^-^ 



(2) 



(3) 



where a dot indicates derivation w.r.t. conformal time, H 
is the Hubble parameter in conformal time, and we have 
defined Jr = df/dR. 

Eq. (3) can be interpreted as an equation of motion 
for a scalar field fu, dubbed scalaron in [21], with an 
effective potential that depends on the density of matter 
and gives the scalaron an effective mass 



l + /j 
Irr 



- R 



l + fR 
"ifRR 



(4) 



There are certain conditions that need to be imposed 
on /(i?) theories in order to avoid instabilities and re- 
produce the correct high curvature regime. Namely, one 
needs fRR > and fRRR ^ 1 to have a stable high 
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today. In a given f{R) model, Bq is in a one-to-one cor- 
respondence with /r - sl heavier scalaron (smaller Bq) 
corresponds to a smaller value of I/rI- In models with a 
LCDM expansion history and — 0.27, they are ap- 
proximately related as Bq « ^Q/r over the range of in- 
teresting values of Generally, models satisfying the vi- 
ability conditions and the local tests of gravity discussed 
above, predict an expansion history practically indistin- 
guishable from LCDM [17-19, 27-29]. Still, the dynamics 
of perturbations can differ significantly [20, 30]. 

To describe the growth of structure and evolution of 
the CMB anisotropics, we expand the f{R) field equa- 
tions to first order in perturbations. The linearly per- 
turbed Einstein equations contain terms in S/r and its 
first and second time derivatives, making them fourth or- 
der in derivatives of the metric (as opposed to the second 
order of GR equations). We will not reproduce the full 
set of exact equations here, as only the Poisson and the 
anisotropy equations will be relevant for our discussion. 
In what follows, we will compare the latter to the equiva- 
lent equations in GR to highlight the modifications intro- 
duced by f{R) theories. We consider scalar perturbations 
in the conformal Newtonian gauge, with Sqqq = — 2a^^ 
and 5gij — ~2a^^Sij, and present all our equations in 
Fourier space. 

Let us start with the anisotropy equation, which is the 
spatial off-diagonal component of the linearized Einstein 
equation. In GR, it reads 
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where a is the anisotropic stress of matter. In /(i?) it 
becomes 
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F 



3a2 (p + P) 
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(7) 



where F = 1 + /i^,. 

The Poisson equation is obtained with a suitable com- 
bination of the time-time and time-space components of 
the Einstein equations that leaves the comoving density 
contrast as the contribution from matter and, in GR, the 
Laplacian acting on the curvature potential as the only 
geometric term. It is, however, common to work with the 
analogous equation for the Newtonian potential, since 
the latter is more directly related to observables given 
that it governs the motion of non-relativistic particles. In 
GR, combining the time-time, time-space and anisotropy 
equations (neglecting the anisotropic stress) , one obtains 
the following algebraic relation: 



In the quasi-static sub-horizon limit, when the terms in 
the square bracket become negligible, Eq. (11) reduces 
to a trivial generalization of (10), i.e. a rescaling of the 
Newton constant by F . 

Overall, the fourth order nature of the linearly per- 
turbed equations in /(i?) theories implies a richer dy- 
namics and potential instabilities, just like at the back- 
ground level. In particular, at early times, i.e. high cur- 
vature and high scalaron mass, the equation for bj^. has 
highly oscillatory solutions [19, 20, 30] with a frequency 
proportional to the scalaron mass. As long as ]b.r > 
(assuming F > 0), these oscillations remain small in am- 
plitude and decay in time (see Sec. II A) ensuring that 
instabilities do not form in the evolution of linear per- 
turbations. 



Linear growth in f{R) and its approximate 
par ametrizat ions 
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where A = S + 31-1/ kv is the comoving density contrast. 
The latter is valid on all scales, as long as the anisotropic 
stress is negligible, and is commonly employed to study 
the clustering of matter. This is the equation that is of- 
ten parametrized along with the anisotropy equation. In 
/(i?), due to the higher order nature of the theory, the 
same combination of Einstein equations leads to a 'modi- 
fied Poisson' equation for which is now dynamical and 
reads: 
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Using Sffi = fuRSR, expanding 5R in terms of the metric 
potentials, and applying the quasi-static approximation, 
Eq. (9) reduces to a simple generalization of (8), which 
will be discussed in Sec. II A. 

It is often preferable to work with a Poisson equa- 
tion for the lensing potential, governing the motion of 
relativistic particles. In GR, when neglecting anisotropic 
stress, the latter reads: 



fc2($ + *) = - 



M 



(10) 



As discussed above, Einstein's equations in /(i?) the- 
ories are quite different from those in GR and can result 
in significantly different dynamics of linear perturbations. 
But modifications of linear growth can also occur in mod- 
els of exotic dark energy and dark matter that are oth- 
erwise based on GR. Hence, we use "MG" to denote not 
only modified gravity, but also, more generally, modified 
growth. 

An arbitrary modification to the dynamics of scalar 
perturbations on linear scales can be encoded into two 
time- and scale-dependent functions, p,{a,k) and ^(a,k), 
generalizing the Poisson (8) and anisotropy (6) equations 
of GR to the following: 
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- =7(a,fc). 



(12) 
(13) 



where any deviation of p, and 7 from unity signals a de- 
parture from the LCDM growth at late times^. The func- 
tions p and 7 do not necessarily have a simple form in spe- 
cific models of MG. Strictly speaking, they parametrize 
solutions of the equations of motions and depend on the 
choice of the initial conditions"^. Nevertheless, they pro- 
vide a consistent framework for searching for departures 
from LCDM. In some theories, including f{R), they can 
assume a simple form on sub-horizon scales. However, 



This equation is again valid on all scales and is conve- 
nient to use when one wants to study weak lensing and 
the Integrated Sachs- Wolfe (ISW) effect. The equivalent 
equation in f{R) is again dynamical: 
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Several other choices of the two functions parametrizing per- 
turbed MG equations, equivalent to fi and 7, can be found in 
the literature. For a summary see [52, 53]. 

This is true also for the conceptually similar but technically dif- 
ferent approach proposed in [54]. Alternative parametrization 
schemes that are more directly related to particular theories have 
also been developed, e. g. in [55-57]. They are either significantly 
more complex, as in [55, 56], or apply to a more limited range of 
models [57]. 
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one can also study them in a model-independent way. 
For instance, in [59-61], a principal component analysis 
of these function was used to forecast and analyze the 
ability of experiments like DES and LSST to constrain 
MG. ' 

Depending on the observables that one considers, it is 
sometimes convenient to parametrize the Poisson equa- 
tion for the lensing potential, rather than the one for 
^' [60, 62, 63]. Namely, one can introduce a function 
E(a, k) via 



($ -t- 4') 



S(a, fc)pA 



(14) 



which is more directly probed by statistical quantities 
derived from weak lensing of distant galaxies. 

In f{R), the Compton wavelength associated with the 
mass of the scalaron (4) sets the range of the fifth force 
interaction, which separates two regimes of dynamics. On 
scales larger than the Compton length, modifications are 
negligible and the dynamics is very close to that in GR. 
Below the Compton scale the growth is enhanced and 
the two metric potentials are no longer equal. This scale- 
dependent behavior can be easily seen from the Pois- 
son (9) and anisotropy (7) equations in the quasi-static 
sub-horizon limit, as discussed at length in [30]. Assum- 
ing that on sub-horizon scales {k ^ aH), the time vari- 
ation of the gravitational potentials is slow compared to 
their variation in space, one obtains a simplified form of 
the equations that can be reproduced substituting the 
following functions into (12) and (13): 



^«(a,fc) 
jQ{a,k) 



1 1 + (4/3)Q 



F 
H 



1 + g 

(2/3)Q 



1 + (4/3)Q 



(15) 
(16) 



where the dimensionless parameter Q [30] is the squared 
ratio of the scalaron Compton wavelength, Ac = 27r/m^^ 
(with nifj^ given by (4)), to the physical wavelength as- 
sociated with k: 
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Ac 
A 



(17) 



Note that the 1/F factor in (15), which enhances the 
growth in a scale-independent way, is small for values of 
l/iji ^ l, favored by local tests of gravity, but is not 
entirely negligible for larger values. Also, for f{R) in the 
quasi-static sub-horizon limit, S(a, k) reduces to a simple 
scale-independent expression, i.e. 



sQ(a) = 1/F 



(18) 



As evident from Eq. (11), there are extra terms that 
would modify this simple relation on larger scales. 

The quasi-static expressions for fi (15) and 7 (16), mo- 
tivated the following general parametrization introduced 



by Bertschinger and Zukin (BZ) in [44] 
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'i,/32,s) and describes 



It contains 5 parameters (ai,a2, 
a transition of both functions from one constant value to 
another, analogously to what happens in the quasi static 
limit of f{R). In fact, the correspondence with Eqs. (15)- 
(16) would be exact if the scalaron Compton length had 
a power law dependence on the scale factor, a. We will 
adopt this as an approximation for now and examine its 
implications later. As discussed earlier, for a given expan- 
sion history, /(i?) theories have a single free parameter, 
fp, or equivalently, Bq. It is easy to express analytically 
cti, a2, Pi and /32 in terms of Bq and one is then left 
with the time dependence, parametrized by s, which can 
be determined numerically once the f{R) model has been 
reconstructed. 

In [45], the authors suggested a slight generalization 
of (19)-(20) to include the pre-factor 1/F whose effect 
can be important for ISW at large Bq. Combining these 
facts, we consider the following expressions, which we will 
refer to as the BZ parametrization: 



A^^^(a,fc) = 



7""(a,fc) 



1 



1 - B^a^-^/Q 
1 + (l/3)SoPa 



1 + {2/3)Bok^a'' 
1 + (l/2)Bofc2a'' 



,(21) 



(22) 



1 + {2/3)Bok^a^ 

where we have introduced a dimensionless wavenumber 
k = k/Ho, such that k = 2997.9 fc Mpc/h. Compar- 
ing (21)-(22) with (15)-(16), it is easy to see that we have 
effectively set A^ « 2'it'^/Hq Bq a^+^. We have also used 
the numerically found approximate relation Bq ss 
in the 1/F pre-factor. 

In principle, the time dependence of the Compton 
length depends on Bq and, strictly speaking, s is not 
an independent parameter. Indeed, the validity of the 
power law assumption was questioned in [58]. On the 
other hand, after it was suggested in [46] that f{R) mod- 
els which closely mimic LCDM correspond to s « 4, 
many authors have used (21)-(22) with s = 4 to place 
constraints on /(i?) theories in terms of bounds on Bq. 
One of the aims of this paper is to investigate the ac- 
curacy of s = 4 and we show in Sec. Ill that it suffices 
to simply fix the value of s to a constant value for the 
range of f{R) models that will be probed with upcoming 
surveys. 

In Sec III, we numerically check the accuracy of the 
quasi-static (/i'^,7'5) parametrization for a range of ffj 
(Bq). We also check how well the time dependence of Q is 
described by a power law in a. Additionally, we estimate 
the theoretical uncertainty in if one were to use ijP^ 
instead of the exact solutions. Before proceeding with 
these tasks, we first revisit analytically the quasi-static 
approximation and the potential role of near-horizon ef- 
fects. 
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B. Quasi-static approximation and near-horizon 
effects 



To obtain the functions and 7 as defined in (15) and 
(16) one applies the quasi-static approximation which 
amounts to neglecting time derivatives of the perturbed 
quantities. Such an approximation is certainly very good 
on scales that are well inside the horizon, for which 
k/aH <C 1, and still in the linear regime; it is less clear 
whether this approximation is safe on near horizon scales. 
It is convenient to work with Sfn and <&+ = ($ + 5')/2 as 
the two geometric scalar degrees of freedom, as in [30]. 

If we perturb the trace equation to linear order, take 
the super-horizon limit and assume a heavy scalaron 
mass, we obtain a damped inhomogeneous harmonic 
equation for S/r, similar to the one discussed in [19]: 



a2 /3(5F 
3Mf \~6i 



l]pS. (23) 



As long as Jrh > 0, the solutions of this equation are 
damped oscillations with frequency and amplitude pro- 
portional to the mass of the scalaron. At early times the 
scalaron is heavy and the oscillations have very high fre- 
quency. As time evolves, the scalaron becomes lighter and 
the frequency decreases, while the amplitude increases. 
Eventually, as the scalaron mass decreases, the damp- 
ing term and the other terms that were neglected in the 
heavy scalaron limit of Eq. (23) become important, and 
the oscillations die off. 

No oscillations are present in the lensing potential. 
Taking the heavy scalaron limit of the equation for $_(-, 
one can see that it decouples from the equation for Sfu, 
becoming: 



2Ml 



kF 



(24) 



and $+ follows closely the evolution it would have in 
GR, except for the \/F factor rescaling the Newton's 
constant. 

In terms of Sfn and $4-, the metric potentials are given 
by $ = $+ + 5fR/2F and * = $+ - Sfii/2F. There- 
fore, they display an oscillatory mode on top of the non- 
oscillatory term $_|_. Their oscillations are out of phase, 
with the same frequency as 6ffj and roughly half the am- 
plitude (since F « 1). 

Oscillations are observed also in the comoving mat- 
ter density contrast, on top of the underlying growing 
mode. They have the same frequency as those in Sfji 
and a relative amplitude which is slightly bigger than 
that of 6fB. This can be easily understood by looking at 
the f{R) Poisson equation (9). In GR, we simply have 
A = -2/cV(3a^^f^)*, vahd on all scales, while Eq. (9) 
contains extra dynamical terms, which are important on 
super-horizon scales; in particular the term c>c T-L^Sfn in- 
duces oscillations on top of the growing mode with an 
amplitude of sa H'^ / {HlQ,jn{a))5 Jr. 



Finally, because of their definition, /i and 7, also dis- 
play oscillations at low k (see e.g. Fig. 2). While the am- 
plitude of the oscillations in 7 increases in time (as it does 
for (5/ij, $, 5* and A), the amplitude of the oscillations in 
/i decreases because A is a growing function. 

These oscillatory modes, generated by dynamical 
terms that are important on super-horizon scales, pose 
certain challenges when one evolves the equations of mo- 
tion numerically. On the other hand, they have practi- 
cally no impact on observables. Namely, they produce 
effects of order H^/mj^, which could be non negligible 
only for larger values of |/^|. Moreover, they only matter 
on the largest scales, and the only observable that could 
show an imprint of these oscillations would be the ISW 
effect. However, the latter is sourced by the time evolu- 
tion of the lensing potential which, as we saw, does not 
oscillate. Therefore, when comparing f{R) theories with 
data, it is safe to work with the {pP ,^'^) parametriza- 
tion obtained in the quasi-static limit, and neglect the 
near-horizon terms. 



III. PARAMETRIZED VS EXACT: 
NUMERICAL COMPARISON 



In this Section, we compare the functions p and 7 ob- 
tained by solving the exact linearized f{R) equations of 
motion to their expressions in the quasi-static approxima- 
tion, (15)-(16), and to the BZ representation, (21)-(22). 
To obtain the "Exact" functions fj^^ and 7^'^, we nu- 
merically solve the full system of linearized Einstein and 
energy-momentum equations, given in [30], to calculate 
A(a, /c), $(a. A:) and ^'(a, fc). Then, we take their ratios, 
as prescribed by Eqs. (12) and (13), to find /Lt and 7. 

Fig. 1 shows contour plots of pF'^ in the [z, k) plane, for 
several values of the parameter We opt to present the 
time dependence in terms of the redshift z, and within 
in the (z, k) domain roughly corresponding to the range 
of linear scales that can be probed by future surveys. 
The main features of jjF''^ are well-described by Eq. (15). 
On scales larger than the scalaron Compton wavelength, 
when Q ^ 1, we have fi = 1 because the growth of 
perturbations is the same as in GR. On smaller scales, 
Q ^ 1, the growth is enhanced by the fifth force and 
/i = 4/3. In addition to this scale-dependent transition, 
the strength of the gravitational interaction is enhanced 
by an overall scale-independent factor 1/F = (l-|-/ij)~^, 
which is very close to unity except for larger values of 
Iffj l - For example, the 1/F enhancement is clearly visible 
in the case of f'^ = — 10~^. 

The behavior of 7^^^ is captured by Eq. (16) and is 
similar to that of p,^^. We do not show plots of j^'^{z, k), 
noting instead that its main difference from /x is the ab- 
sence of the 1/F factor, while the transition is from 7 = 1 
on large scales to 7 = 1/2 below the Compton wavelength 
scale. 
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fi and 7*^ 



In the quasi-static regime of linear sub-horizon per- 
turbations, the f{R) modifications to the Poisson and 
anisotropy equations take the form of fjfi (15) and 
7*^ (16). The background functions fn and Jrh appear- 
ing in the definition of Q can be found numerically by 
solving the f{R) Friedmann equation (2) for a given ex- 
pansion history. 

To investigate the accuracy of the quasi-static approx- 
imation, we estimate the relative difference between jjfi 
and jj^^ for values of from 10^^ to 10^® and show the 
results for some of these models in Fig. 2. We find that the 
differences decrease as \f^\ gets smaller. For < 10^^ 
the average difference is about 0.01%, while the maxi- 
mum difference is about 2.0%. The largest errors occur 
for f'j^ = — 10~^, where the maximum difference is about 
7%. 

As the plots in Fig. 2 show, the main difference is due 
the oscillations in time present on large scales in jj^^ and 
that are not present in /i'^. As discussed in Sec II B, 5 fa 
has an oscillatory behavior that translates into oscilla- 
tions in A, <I> and ^I', but not in $ + ^I*. The oscillations 
are most prominent in the homogenous {k — 0) limit and 
die off on smaller scales. While they dominate the rela- 
tive difference between /z*^ and they are actually 
quite small and, as discussed earlier, they do not have 
any observable consequences. 

In addition to the oscillations, the /]j = — 10~^ case 
shows visible differences at late times. These are due to 
the fact that some of the neglected terms in the perturbed 
equations are multiplied by time derivatives of fu and 
/rRj which can become non- negligible for larger values of 
|/^|. Such high values of |/^| however, are already ruled 
by the existing cluster abundance data [64, 65] which, 
combined with other cosmological probes, places a bound 

of \m < 10-3. 

The examination of the differences between ^'^ and 7^^^ 
yields very similar conclusions. Overall, we see that the 
quasi-static expressions (15) and 7*^ (16) are an ex- 
cellent representation of f{R) solutions for < IQ-^. 



B. 



and 7^^ 



The BZ parametrization of Eqs. (21) and (22) is in- 
spired by the quasi-static approximation corresponding 
to and 7*^, but it contains one more degree of ap- 
proximation - it assumes a particular form for the time 
dependence of Q. Namely, it sets 



6^0 Irr 



(25) 



In [46] it was suggested that s « 4 for a wide range of 
models. We want to test the accuracy of this assumption 
and the impact the choice of s makes on /(i?) constraints 
derived using BZ. 



rO 
JR 


-10-1 


-10-^ 


-10-^^ 


-10-" 


-lo-'^ 


-10-"^ 


S 


3.6 


3.6 


3.6 


3.5 


3.3 


3.2 



TABLE I. The best fit values of s, obtained by fitting a" to 
the exact time dependence of Q using Eq. (25), for six values 
of/S- 



We determine numerically the time-dependence of Q, 
and fit a'^ to the expression in Eq. (25). Table I shows 
results obtained by fitting over the redshift range < 
z < 3 for models with 10-'^ < I/r] < 10"^ We see that 
s = 4 is a reasonable choice for \f^\ > 10^^, but for 
smaller \ffi\, values closer to s = 3 are preferred. Thus, a 
fixed s, and s = 4 in particular, does not appear to hold 
for all However, while there is a clear dependence of s 
on f'^, we will argue shortly that such differences should 
not be taken very seriously. 

As already pointed out in [58] , the time dependence of 
Q is not exactly a simple power law. This is because the 
onset of the accelerated expansion at late times changes 
the evolution of the scalaron mass compared to that in 
the matter domination epoch. Fig. 3 shows the numeri- 
cally reconstructed time dependence of Q for f'^ = —10-^ 
and —10-^, along with plots of a* with s = 3, 3.5 and 
4. Clearly, s = 4 is a better fit at early times, during 
the matter dominated epoch, while there is a transition 
towards s = 3 near the matter-A equality. This transi- 
tion can be modeled using a phenomenological function 
of the form s = 3.5 — 0.5tanh[(a — 0.45)/0.3], in which s 
transits from 4 to 3 around the time of the matter to A 
transition. It turns out, however, that whether one uses 
s = 4, s = 3.5, or the phenomenological function, the 
differences between fi^^ and /i^^ remain roughly of the 
same order. This can be seen comparing Figs. 4, 5 and 6, 
where we show the percent difference between fi^^ with 
s ^ 4 , s = 3.5 and s = 3.5 - 0.5 tanh[(a - 0.45) /0.3] and 
the numerically found 11^^. One can see that, while the 
differences are slightly larger than those between /^"^ and 
fi^^, they follow the same general trend. The most visi- 
ble difference is due to oscillations in at small k. As 
expected, using the phenomenological functions results 
in the smallest discrepancies. However, in all three cases, 
the differences are less than 2% - 3% for < 10-^. 
The same conclusions hold for 7^^. 

To be more quantitative, we use a Fisher matrix for- 
malism to estimate the variance in Bq given the theoret- 
ical error associated with using //^^. The Fisher matrix 
element corresponding to Bo is calculated as 



BqBq 



BZ 



dBo 



(26) 



where Cij is the covariance matrix for /i, and /if ^ denotes 
the value of /i^^ at a particular (z, k) pixel on our grid. 
We pixelate the (z, k) domain into 20 bins in z and 20 
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FIG. 3. Time dependence of Q for — — 10~^ (orange dotted-dashed line) and = —10"'* (purple dashed line), compared 
to the power laws a* (blue crossed line), a^'^ (red dotted line) and a!^ (green double dotted-dashed line) as well as to the 
phenomenological function s = [3.5 — 0.5 tanh(a — 0.45)/0.3] (indicated as 'tanh fit', in solid black). The phenomenological 
function can follow better the transition from a* at early times to a?' at late times. Note that Bo ~ — 6/^.. 
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bins in k, using the same binning configuration as the 
one used for the PCA analysis in [59, 60]. The partial 
derivatives are calculated analytically using Eq. (21). We 
take the covariance matrix to be diagonal, 

C,, = (/.r-Mr)<5,, , (27) 

which effectively treats fif^ as a biased "measurement" 
of /if ^ that is made independently in each bin. We then 

take ^ ^BoBo ^ rough estimate of the variance in Bq. 
The variances computed from Eq. (26) are two orders 
of magnitude smaller than their corresponding fducial 
values of Bq, i.e. 



(28) 



for 10^^ < Bq < 0.1. The forecasted experimental con- 
straints on the /i (and 7) from LSST [60], based on the 
same pixelation scheme, are much larger than the system- 
atic theoretical errors due to using /i^^ and 7^^. Hence, 
we conclude that the BZ approximation can be safely 
used to put constraints on the f{R) models. 
We also checked that there is essentially no significant 
difference between using s = 4 versus s = 3.5, or the 
phenomenological fit function in the computed values for 
the variance of Bq. This insensitivity of the uncertainty 
in Bq to the value of s is because the dependence of fi 
and 7 on Bq is much stronger than its dependence on 
the time variation. This can be seen by using Eq. (21) to 
obtain 



d^/ds 
W/dBQ 



_Bo In a , 



(29) 



which is small for Bq < 0.1 and 0.1 < a < 1. Similar 
conclusion hold for 7. 

As an additional check, we have implemented and 
7^^^ in MGCAMB [47] for several representative values 
of 10~^ ^ \ffi\ ^ 10~^ and compared the output for the 
CMB and linear matter power spectra to that obtained 
using /i^^. We found that, as expected, the differences 
are too small to be of relevance. 



IV. SUMMARY 

We have examined in some detail the parametrizations 
commonly employed to evolve linear scalar perturbations 
in f{R) models and place cosmological constraints. In 
particular, we have carefully examined the validity of the 



quasi-static approximation and of the parametrization in- 
troduced in [44], (BZ), and adapted to /(i?) in [45]. 

After reviewing the main features of the linearized Ein- 
stein equations in f{R), we have analyzed near horizon 
effects and provided an analytical insight on the oscil- 
lations that are observed at early times on large scales. 
Upon showing that the oscillations originate from the 
scalaron, and linking their amplitude and frequency to 
the mass of the latter, we have argued that they have 
no observable signatures. We have then numerically an- 
alyzed the accuracy of the quasi-static approximation, 
finding that on linear sub-horizon scales the functions 
/i*^ and 7*^ describe very closely the modified Poisson 
and anisotropy equations. 

Furthermore, we have analyzed the BZ parametriza- 
tion [44], in the form adapted to f{R), and in particular 
the validity of its approximation of the time-dependence 
of Q as a power law a". We have found that there is not 
a fixed value of s valid for all designer f{R) models, but 
rather the best fit value of s depends on f'^ and the range 
of redshifts considered. A comparison with the exact time 
dependence of Q shows that s w 4 is a good fit for early 
times, while s ~ 3 fits better at late times. The transition 
of the value of s can be modeled using a phenomenolog- 
ical function that takes into account the transition from 
matter to A dominated era. However, we have shown that 
fi and 7 are relatively insensitive to the variation of s so 
that s — 3.5, s = 4, or the phenomenological fit lead to 
a theoretical systematic error in fi (and 7) that is much 
smaller than the observational uncertainty from a future 
survey like LSST. Hence, the BZ parametrization can be 
safely used for deriving constraints on /(i?) models from 
upcoming large scale surveys. 
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